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ABSTRACT 


In  this  paper  a simple  inverse  filter  is  developed  to  suppress  multiple 
reflections  from  a normal  incidence  synthetic  seismogram.  The  filter  is  de- 
veloped by  means  of  Mendel's  Bremmer  Series  Decomposition  (Refs.  1 and  7)  - 
the  state  space  model  for  the  complete  response  is  decomposed  into  primaries, 
secondaries,  tertiaries, . . . , etc.,  models  which  generate  only  primary,  se- 
condary, tertiary,  ...»  etc.  multiple  reflections,  respectively.  The  equation 
formulations  are  based  on  the  operator  description  of  state  space  model  of 
the  layered  media  (ref.  2). 

This  filter  removes  all  the  multiple  reflections  from  the  seismogram 
which  are  ever  reflected  off  of  the  surface,  inside  the  layered  media,  and 
is  especially  effective  when  the  surface  reflection  coefficient  is  rela- 
tively large  as  in  most  geophysical  situations.  A recursive  scheme  of 
applying  this  filter  to  the  consecutive  subsystems  of  the  layered  media  is 
also  developed.  It  generates  successive  approximations  of  the  primaries  of 
the  system,  the  final  result  being  the  pure  primaries. 


I.  Introduction 

In  seismic  profiling  of  layered  media,  primary  reflections  - compres- 
slonal  waves  reflected  directly  from  the  interfaces  of  the  layers  - contain 
the  ultimate  information  for  the  determination  of  the  subsurface  structure 
of  the  media.  Minimization  of  the  contribution  from  other  undesirable  events, 
such  as  noise,  surface  waves  and  ghost  reflections,  to  maximize  the  contri- 
bution from  primary  reflection  has  been  one  of  the  fundamental  problems  in 
the  analysis  of  seismic  traces.  Another  and  particularly  troublesome  type 
of  Interference  is  that  from  multiple  reflections.  There  are  some  special 
classes  of  multiple  reflections  which  exhibit  significant  amplitudes  so 
that  they  occur  well  into  the  seismogram  and  affect  its  appearance.  They 
are  most  likely  to  happen  when  the  layered  system  involves  one  or  more 
strong  reflectors.  They  occur  through  reinforcement  of  several  multiple  re- 
flections. A general  theory  has  been  developed  to  locate  and  identify  various 
kinds  of  multiple  reflections  and  reinforced  events  between  them  (ref.  6). 

In  the  analysis  of  multiple  reflections,  which  represent  a highly  complicated 
mechanism  of  layered  media,  Mendel's  Bremmer  Series  Decomposition  (refs.  1 
and  7)  has  been  found  to  be  useful,  which  makes  it  possible  to  study  each 
of  the  decomposed  multiple  reflections  such  as  primaries,  secondaries,  ter- 
tiaries,  ...,  etc.,  separately.  The  state  space  models  for  lossless  waves 
in  layered  media  which  are  described  by  the  wave  equations  and  boundary 
conditions  have  been  developed  by  Mendel,  et.  al.  (refs.  2 and  3).  The  models 
are  for  non-equal  one-way  travel  times  and,  by  the  nature  of  non-uniformity, 
represent  a special  class  of  equations  with  multiple  delays  which  are  referred 
to  as  causal  functional  equations.  Based  on  these  equations  Mendel  (ref.  1) 
has  proven  the  truth  of  the  following  decomposition  of  the  solution  to  the 
lossless  wave  equation  in  layered  media  : the  complete  output  from  a K-layer 
media  system,  which  is  comprised  of  the  superposition  of  primaries,  secon- 
daries, tertiarles,  etc.,  can  be  obtained  from  a single  state  space  model 
of  order  2K  - the  complete  model  - or  from  an  infinite  number  of  models,  each 
of  order  2K,  the  output  of  the  first  of  which  is  just  the  primaries,  the  out- 
put of  the  second  of  which  is  just  the  secondaries,  etc.  This  decomposition 
of  the  solution  to  the  lossless  wave  equation  into  physically  meaningful  con- 
stituents (i.e.,  primaries,  secondaries,  etc.)  is  called  a Bremmer  Series 
Decomposition,  after  Bremmer,  who  in  1951  (ref.  4)  established  a similar 
decomposition. 
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In  Chis  paper,  by  means  of  Bremmer  Series  Decomposition,  we  develop 
some  methods  to  approximate  the  primary  reflections  of  a layered  system 
from  the  normal  incident  synthetic  seismogram  generated  by  the  system  for 
the  given  information,  such  as  surface  reflection  coefficient  and  input 
waveform,  etc.  From  these  a simple  inverse  filter  is  derived  which  suppresses 
multiple  reflections.  The  filter  removes  all  the  multiple  reflections  from 
the  seismogram  which  are  ever  reflected  off  of  the  surface  inside  the  media. 
This  filter  is  especially  effective  when  the  surface  reflection  coefficient 
is  relatively  large  as  in  most  geophysical  situations.  A recursive  method 
is  also  developed  which  applies  the  filter  to  the  consecutive  subsystems 
of  the  layered  media  from  the  top  to  the  bottom  layer,  generating  succes- 
sive approximations  of  the  primaries  of  the  system. 

In  the  first  few  sections  in  this  paper  (from  section  II  through 
section  V),  we  review  the  state  space  model,  the  Bremmer  Series  Decompo- 
sition and  its  operational  and  transfer  functional  descriptions.  The  inverse 
filter,  which  is  denoted  by  throughout  this  paper,  is  developed  in 
Section  VI.  In  Section  VII,  another  inverse  filter  is  introduced  which  is 
derived  from  a further  approximation  of  primaries  with  some  additional 
information  of  the  layered  system.  The  effects  of  those  filters  and  their 
relationship  are  examined  in  Section  VIII,  from  which  the  recursive 
method  of  applying  filter  F^  to  the  subsystems  of  the  media  is  developed 
in  Section  IX.  For  illustrative  purpose,  a three  layer  example  is  used 
throughout  the  analysis. 
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II.  A State  Space  Model 


A state  space  model  for  the  system  of  K-layered  media,  depicted  in 
Figure  1,  is  derived  under  the  following  modeling  assumptions;  (1)  plane 
compressional  waves,  (2)  horizontally  stratified  nonabsorptive  layers  of 
different  travel  times,  and  (3)  normally  incident  waves.  Each  layer  is 
characterized  by  its  one  way  travel  time,  t^,  velocity  and  normal  in- 
cidence reflection  coefficient  ^ (i  - 1,2,. . . ,K).  Additionally,  interface-0 
denotes  the  surface  and  is  characterized  by  reflection  coefficient  rQ.  We 
adopt  the  convention  of  calling  the  layer  below  layer  K the  basement.  In 
Figure  1,  m(t)  and  y(t)  denote  the  input  (e.g.,  seismic  source  signature 
from  dynamite,  alrgun,  etc.)  to  the  layered  earth  system  which  is  applied 
at  interface-0,  and  the  output  (i.e.,  ideal  seismogram)  of  the  system 
which  is  observed  at  the  surface  respectively. 

The  compressional  waves  within  the  k-th  layer  are  identified  by  two 
states  u^  and  d^,  which  denote  the  upgoing  and  downgoing  waves  in  the 
k-th  layers,  respectively.  At  present  time  t,  uk  is  defined  at  the  top 
of  layer  k,  whereas  d^  is  defined  at  the  bottom  of  layer  k,  as  shown  in 
Figure  2.  To  develop  the  state  equation  model  we  direct  our  attention  at 
the  intersection  point  of  the  ray  diagram  and  apply  superposition  to 
obtain  the  following  equations  for  signals  uk(t  + Tk>  and  d^^  (t  + Tk+^) , 


which  leave  that  point; 

“ rkdk(t>  + ^-rk>  uk+i(t)  (la) 

dk+l^t  + Tk+l^  * (1  + rk^  dk(t)  " rk  uk+l(t)  (lb) 

These  equations  are  applicable  for  k-1,2 , . . . ,K-1.  At  the  surface 
(Figure  3a),  we  obtain 

dl(t  + Tl)  " ~r0  ul(t)  + <1  + r0>  m(t)  (2) 

y(t)  - (i-r0)  Uj(t)  (3) 
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and  , at  the  K-th  interface,  we  assume  that  uR+^(t)*0,  to  obtain 
(Figure  3b) 

uR(t  + tr)  »*  rR  dR(t)  (4) 

Signal  y(t)  in  Eq.  (3)  is  the  measurable  system  output.  Signal 
dR+i (t) , which  is  the  downgoing  wave  in  the  basement,  is  also  a system 
output;  but,  since  it  cannot  be  measured,  we  shall  ignore  it  in  the 
following  analyses. 

It  is  convenient  to  group  Eqs.  (la),  (lb),  (2)  and  (4)  in  a layer 
ordering,  as  follows  : 

dl(t  + Ti)  - -rQ  u^t)  + (1  + rQ)  m(t) 

ul^t  + TP  “ rl  + u2  Ct) 

V'+V  ' (1  + tJ-l)  Vl(°  ' fJ-l  "J(t)j  j.2,3 K-l 

Uj^t+Tj^  ” ^ + “j+l^ 

dK(t  + TK^  * (1  + rK-P  dK-l(t)  ' rK-l  UK(t) 

uK(t  +tk)  - rR  dR(t)  (5) 


Equations  (5)  and  (3)  together  represent  the  state  equation  model  for 
the  complete  output  y(t);  hence,  they  are  referred  to  in  the  sequel  as 
the  complete  model. 

Comment . Equations  (5)  and  (3)  can  be  expressed  in  more  compact  no- 
tation by  introducing  the  following  2K x 2K  matrix  operators  : 


- diag  (z1,z1,z2,z2,...,zk,zk). 


(6) 


where  z^  is  a scalar  operator  used  to  denote  a sec.  time  delay 
(i.e.,  z^t)  - f(t-Ti)).  Let 

jC( t)  - col  (u1(t),d1(t),u2(t),d2(t),...,uK(t),dR(t))  ; 

then  Eqs.  (5)  and  (3)  can  be  written,  as 

SLl J2C(t)  - A ^(t)  + b m(t)  (7a) 

y(t)  - tf^(t)  (7b) 


* In  refs.  1,2, 3, 6 and  7,  the  output  equation  was  defined  as 
y(t)  ■ (l~ro)  u (t)  + rQ  m(t)  including  the  direct  reflection 
term  r m(t;.  Here  we  neglect  the  direct  reflection  term. 
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where  the  explicit  structures  of  A,  b and  c^  can  be  deduced  directly  from 
the  former  equations.  State  Eq.  (7a)  is  a dynamical  equation  with  multiple 
time  delays.  It  is  not  a differential  equation,  nor  is  it  a finite-dif- 
ference equation.  We  shall  refer  to  it  as  a causal  functional  equation. 

It  is  linear  and  time- invariant,  and,  as  is  the  case  with  delay-time 
systems,  requires  initial  value  information  over  initial  intervals  of  time. 
Consider  the  k-th  layer  (Figure  2).  Then  dfc(t)  is  equal  to  zero  until 
t»T^  + t ^ + ...+1^,  and  u^  is  equal  to  zero  until  t * + ^ +•••+  2t^. 

These  facts  are  true  for  all  k«l,2,...,K;  i.e., 

r JL 

d,(t)  - 0 V t e [0  , £ t ) (8a) 

J i-1  1 

and 

r J 

Uj(t)  - 0 V c E 1°  * Ti  + Tj)  (8b) 

where  j • 1,2, . . . ,K. 

For  more  detailed  discussions  about  the  derivation  and  structure 
of  the  state  space  model,  the  reader  is  referred  to  refs.  2 and  3. 
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III.  A Bremmer  Series  Decomposition 


Since  the  complete  output  y(t)  is  a superposition  of  primaries, 
secondaries,  tertiaries,  etc.,  it  can  be  written  as 

oo 

y(t)  - £ • o) 

j-1 

where  (t)  denotes  the  j-ary  reflection  components  of  y(t).  In  this 
paper  we  just  summarize  this  decomposition  in  the  following  theorem, 
given  without  proof  (see  refs.  1 and  3 for  the  proof). 

Theorem 

The  complete  output,  y(t),  from  a K-layer  media  system  can  be  obtained 
from  a single  model  of  order  2K  - the  complete  model,  given  by  Eqs.  (5) 
and  (3)  - or  from  an  infinite  number  of  models,  each  of  order  2K,  inter- 
connected as  shown  in  Figure  4.  The  primaries  model,  which  generates  only 
primary  reflections,  and  the  n-aries  models,  each  of  which  generates  only 
n-ary  reflections,  where  n-2,3,...,  are  defined  by  the  following  : 

(a)  Primaries  model 


and 


dj  i/t+Ti)  = (l+rQ) 

Ul,l(t  + Tl)  " rL  dl,l(t)  + (1_rl)  UI,2(t) 

dl,j(t  + V " (1+rj-l}  dl,j-L(t)  ( 
Ul.j(t  + V “ "j  dl,j(t)  + (l“rj)  ul.j+l(t>* 

dl,K(t  + TK)  " (1+rK-I)  dl,K-l(t) 

Ul,K(t  + TK)  ’ rKdI,K(t> 


yi(t)  “ (1-r0)  ui,l(t) 


j - 2 , 3 , . . . , K- 1 


(10) 


(ID 
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(b)  n-aries  model 


and 


dn,I<t  + Tl) 


-ro  Vi,i(t) 


u .(t  + T ) - r.  d ,(t)  + (1-r  ) ,(t) 

n,  l l i n,  l l u , i 

d .(t  + T ) * (1  + r ) d (t)  - r u (t) 

n,j  j j-1  n,j-l  j-1  n-l,j 

u (t  + T ) - r d (t)  + (1-r  ) u ...  (t) 

n,j  j J n,j  j n, j+1 


d„,K<t  + ,lc' 


(1  + tK-l>  dn,K-I(t)  - rK-l  Vl,K(t) 


%,K(t  + V " rK  d»,K(t) 


y„(t> 


(1-r0) 


j - 2,3,...,K-1 


(12) 


(13) 


for  n ■ 2,3, . 


Additionally,  y(t)  is  given  by  Eq.  (9).  In  these 


equations,  d . , u . and  y denote  n-ary  dovmgoing  and  upgoing 
n,j  n,j  n 

states  and  the  n-ary  reflection  portion  of  the  complete  output, 
respectively. 


The  primaries  model  is  obtained  from  the  complete  model  given  by 
Eq.  (5)  by  deleting  the  term  ~rj_j  (t)  in  the  equations  for  d^(t  + Tj), 
j-l,2,...,K.  This  is  done  to  truncate  the  multiple  reflections  higher 
than  secondary  reflections  which  are  due  to  the  upgoing  waves  reflec- 
ting off  of  the  top  of  the  layers.  The  n-aries  model  is  obtained  in  a 
similar  manner,  by  successively  subtracting  all  the  j-aries  models 
(where  j ■ 1,2, . . . ,n-l) , from  the  complete  model  to  obtain  a residual 
model,  and,  by  then  deleting  the  terms  which  truncate  multiple  reflec- 
tions higher  than  n-ary  reflections  in  that  residual  model.  The  reader 
who  is  interested  in  the  detailed  derivation  of  the  n-aries  model  is 
referred  to  ref.  1 or  7. 

Instead  of  a formal  proof  of  the  above  theorem,  we  demonstrate  its 
validation  through  a three-layer  simulation  (Figures  5 and  6).  Figure  5 
depicts  a three  layer  media  which  can  be  associated  with  a bright  spot 
phenomena,  because  of  the  thin  low  velocity  layer  which  is  sandwiched 
in  between  the  two  thick  high  velocity  layers.  We  attribute  no  other 

geological  plausibility  to  this  example.  For  layer  1,  V ■ 7,500  ft/sec 

3 1 3 

and  p^*2,2  gm/cm  ; for  layer  2,  “ 5,500  ft/sec  and  p^m  1*6  gm/cm  ; 

for  layer  3,  V3“Vj  and  p3*pi»  ^or  t'ie  t>asement>  * 12,000  ft/sec 

and  was  approximated  by  the  1 /A  - power  law,  0,23  (ref . 5). 
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Figure  6 depicts  the  complete  response  as  well  as  the  primaries, 
secondaries  and  some  of  the  tertiaries  (through  2 seconds)  which  were 
obtained  via  simulation  of  Eqs.  (10)—  ( 13) . 

In  this  example  we  observe  that  the  superposition  of  the  first 
three  terms  in  the  Bremmer  series  decomposition  is  a good  approximation 
to  the  complete  response.  In  many  geophysical  situations,  where  reflec- 
tion coefficients  are  quite  small,  the  decomposition  can  be  truncated 
after  secondaries  or  tertiaries;  hence,  the  Bremmer  series 
decomposition  also  represents  a way  to  approximate  the  solution  to  the 
wave  equation.  This,  together  with  the  fact  that,  by  means  of  the  de- 
composition, it  is  possible  to  deal  with  each  of  the  multiple  reflec- 
tions separately,  simplifies  complicated  problems  in  the  analyses  of 
a seismogram. 


IV.  Operator  Descriptions  for  Bremmer  Series  Decomposition  and  Layer 
Transition  Matrix  H 


The  state  space  equations  (10)  for  the  primaries  model  can  be  expres 
sed  in  a compact  way  by  introducing  the  following  KXK  matrix  operator, 

„ A 


diag  (z ^ , z2 » • • • , z^) , 

where  z^  is  a scalar  operator  used  to  denote  a sec.  time  delay 
(i.e.  z^  f (t)  ■ f (t  - t^)) , and,  by  reordering  the  equations  in  such  a 
manner  that  all  downgoing  states  are  grouped  together  and  all  upgoing 
states  are  grouped  together.  Let 


(14) 


d^  (t)  - col 

(dl, 1 (t) * dl,2(t)*  dl,K(t)) 

(15a) 

Uj  (t)  ■ col 

(“l , 1 (t) ' ul,2(t) ul,K(t))* 

(15b) 

Then,  Equations  (10)  can  be  written,  in  partitioned  form  as  : 
Z'1  dL(t)  - Al  d^t)  + g m(t) 


,-l 


where 


(1  + rL  ) 
0 
0 


Z u^t)  - A3  dj(t)  + ux(t) 


0 

0 

'(l  + r2) 


0 

0 

0 

(l  + r3) 


0 

0 

0 

0 


0 

0 

0 

0 


...  ( 1 + r ^ ) 0 


(16a) 


(17) 


0 

0 

0 

• 

• 

0 

0 


1 


(1-r,) 

0 

0 

• 

• 

0 

0 


A^  m diag  (r ^ ,r 2 > • • • « r^) 
0 


(18) 


(l-r2) 


0 

0 


0 

0 

(l-r3) 


0 

0 


0 

0 

0 


«-rK_i> 


(19) 


and 
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(20) 


g - col  ( 1 + rQ , 0,  0,  ....  0)  . 


Matrices  Aj,  A3  and  A^  are  K x K,  and  g.  is  K*  1.  Equations  (16)  can  be 
solved  for  u^(t)  as 


v-1 


u^t)  - (I-ZA4)  4 ZA3(I-ZA1)"1  Z*  m(t)  . 

In  addition  to  Eq.  (21),  we  have  the  observation  Eq.  (11),  which  can 
be  written  as 


(21) 


or 


where 


y^(t)  ■ (0,0, •••,0  i (1  - Tq) ,0, . . . ,0) 


yx(t)  = h'  u^t) 


h'  ■ ((1  - Tq) ,0,. . . ,0) , which  is  lx  K . 


d:(t) 


u.  (t) 


(22) 


(23) 


In  the  same  manner,  Eq.  (12)  for  n-aries  model  can  be  expressed 


as 


."I 


Z d^(t)  - A1  d^(t)  + A2  u^^t) 


(24a) 


Z 1 u (t)  - A,  d (t)  + A.  u (t)  (24b) 

n j — n 4 n 

where  A2  is  a K x K matrix  given  by 

Aj  “ diag  (-Tq , — r j , . . . , ~ r ^ ) (25) 

We  notice  that  in  Eq.  (24a),  the  input  to  the  system  of  n-aries  model 

is  the  upgoing  waves,  u ^ ^ (t)  of  the  (n-l)-aries  model.  Solving  Eqs. 

(24a)  and  (24b)  for  u (t),  we  find  that 

— n 

^(t)  - (I-ZA4)-1  ZA3(I-ZA1)‘l  ZA2  u^t)  . (26) 

Additionally  (13)  can  be  expressed  as 

y„(t>  - h’  U (t)  (27) 

n — n 

where  h’  is  given  by  (23) 
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We  observe  that  the  same  matrix  term  (I-ZA^)-1  ZA^  (I  - ZA^)”*Z 
appears  in  Eqs.  (21)  and  (26).  Let 

H - (I  - ZA4)_1  ZA3(I-ZA1)"1  Z . ( 


Then,  Eqs.  (21)  and  (26)  can  be  expressed  in  terms  of  H,  as 

Uj(t)  * H & m(t) 

Un(t)  - H A2  Un_1  (t) 


(29a) 

(29b) 


These  are  the  recursive  equations  for  u^t).  In  non-recursive  form, 
we  can  write  Eq.  (29b)  as 

u^t)  - (HA2)n-1  uL(t)  . i 


Explicit  expressions  for  matrix  H in  terms  of  reflection  and 
transmission  coefficients  with  multiple  delay  operators  were  obtained 
in  ref.  6.  Matrix  H is  KxK,  for  a K-layer  system,  and  is  given  by 


H - {h^j ) 


‘ K I / 

.»?/"  Vl  Vi  pi-l  Vl 


(31a) 


for  i > j,  i,j  * 1,2,. ..,K 


Hij  “ hji 


(31b) 


p q 


where 


1 2 1 2 1 

e - n z , s - n (i-o  . p,  - n 0-0 

ji-i  1 2.- 1 1 2-1  1 


qi  “ Jl(1  + ri)*  and 
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V.  Transfer  Function  Representations  of  Bremmer  Series  Decomposition 


Taking  the  Laplace  Transformation  of  Eqs.  (10),  (11),  (12)  and  (13), 
and  introducing  the  following  K*K  matrix. 


A r -T,s  “T-s 

Z(s)  - diag  j^e  1 , e 2 e J 


we  obtain  equations  of  the  same  structure  as  Eqs.  (16a),  (16b),  (24a) 
and  (24b);  i.e., 

Z_1  (s)  D^s)  - Al  D^s)  + * M(s) 

Z_1(s)  Uj  (s)  - A3  D^s)  + A4  Uj  (s) 


(33) 

(34a) 

(34b) 


for  primaries  model,  and 
,-l. 


Z (s)  D (s)  - A D (s)  + A,  0 (s) 

l ti  2 ~ n-l 

Z_1(s)  U (s)  - A D (s)  + A.  U (s) 

J —n  4 — n 


(35a) 

(35b) 


for  n-aries  model,  where 


D (s)  - col  (D  (s),  D (s) r(s)) 

n,  1 n,Z  n,K 

U (s)  - col  (U  (s),  U (s) U (s)) 

n,  i n,z  n,K. 


n * 1 ) 2 ) • • • ) K. 


(36a) 

(36b) 


It  is  straight  forward  to  show  that  the  Laplace  transforms  of  the 
rest  of  the  equations  in  Section  IV  also  remain  the  same  in  their  expli- 
cit structures.  Especially,  the  Laplace  transforms  of  Eqs.  (29a),  (29b) 
and  (30)  are  given  by 

U^s)  « H(s)  & M(s)  (37a) 

1^(8)  - H(s)  A2  U^s)  (37b) 

and 

1^(8)  - (H(s)  A2)n_1U1(s)  (38) 

respectively,  where 

H(s)  - (I-Z(s)  A4)_1Z(s)  A3(I- Z(s)  A1)_1.Z(s)  (39) 

which  is  the  Laplace  transformed  version  of  the  layer  transition 
matrix  H given  by  Eq.  (28). 
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In  addition,  the  output  equation  is  given  by,  from  Eqs.  (22)  and 


(27), 


n * 1,2, ...  ,K. 


Yn(s)  - h'  UJs) 


Z-transformation  is  also  applicable  in  our  Bremmer  series  decom- 
position if  we  let 

r^-n^T  , i-  1,2,. ...K 

where  T is  the  sampling  interval.  As  we  did  in  Laplace  transfor- 
mation, defining 


Z (z)  * diag  jV**,  z z KJ  , 


we  obtain  z-transformed  version  of  Eqs.  (16a),  (16b),  (24a)  and 
(24b);  i,e. , 

Z_1(z)  D^z)  - Al  Dl(z)  + £ M(z) 

Z_1(z)  U^z)  - A3  Dx(z)  + A4  U^z) 
for  primaries  model,  and 

Z(z)_1  D^z)  - Al  D^z)  + A2  U^z) 

Z(z)"1  U^z)  - A3  D^z)  + A4  U^z) 


for  n-aries  model,  where 


D (z)  - col  (D  .(z),  D ,(z),  ....  D „(z)) 
n,  l n , z n , K. 

U (z)  - col  (U  ,(z),  U ,(z),  ....  U „(z)) 
n,  i n,z  n,K 


(43a) 


(43b) 


(44a) 


(44b) 


(45a) 


(45b) 


N-  1,2,.  ...K. 


The  z-transforms  of  the  rest  of  the  equations  in  Section  IV  also 
remain  the  same  in  their  explicit  structures.  Consequently,  the  operational 
representations  in  Section  IV,  the  Laplace  transforms  and  the  z-transforms 
are  readily  transferrable  to  each  other  merely  by  changing  the  notations 
of  corresponding  variables. 
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The  corresponding  z-transforms  of  Eqs.  (29a),  (29b)  and  (30) 
are  given  by 

U^z)  - H(z)  & M(z)  (46a) 

V2>  - H(z)  A2  Un_1  (z)  (46b) 

and 

Vz)  - (H(z)  A2)n'1U1(z)  (47) 

respectively,  where 

H(z)  - (I-Z(z)A4)_1Z(z)A3  (I-Z(z)A1)'1  Z(z).  (48) 

In  addition,  the  output  equations  are 

Yn(z)  - h'  U^z)  (49) 

n ■ 1,2, .. . ,K. 


It  is  not  difficult  to  show  that 

H(s)  - H ^ ~Tis 
z^-*  e 

H(z)  - H ni 

z ]_  •¥  Z 

i*l,2,...,K.  We  refer  to  H(s)  and 


(50a) 

(50b) 


layer  transfer  matrix. 


VI.  Derivation  of  the  Inverse  Filter  for  Suppressing  Multiple  Reflections 

In  the  following  analysis  we  restrict  ourselves  to  the  system  of 
layered  media  for  which  the  modeling  assumptions  given  in  Section  II 
apply.  We  will  use  the  z-transform  representations  of  our  state  space 
model  in  all  derivations  throughout  this  analysis;  for,  it  is  usually 
convenient  for  simulation  purposes. 

Our  objective  is  : given  a synthetic  seismogram  with  some  infor- 
mation about  the  system,  (such  as  estimated  reflection  coefficients  and 
input  waveform)  to  suppress  multiple  reflections  from  the  seismogram  so 
as  to  maximize  the  contributions  from  the  primaries.  In  our  state  space 
model,  it  is  equivalent  to  say  that  given  Y(z),  the  complete  response 
of  the  layered  media,  we  want  to  find  an  inverse  filter  which  gives 
Y^(z),  the  primaries  (the  output  of  the  primaries  model),  or  an  appro- 
ximation of  Yj(z)  from  Y(z). 

Suppose  the  reflection  coefficient,  of  the  surface  and  the  in- 
put waveform,  M(z),are  both  known,  and  that  the  seismogram  Y(z)  (the 
output  of  the  system  generated  by  M(z))  is  also  known.  We  are  going  to 
suppress  multiple  reflections  from  Y(z)  as  much  as  possible  making  use 
only  of  Tq,  M(z)  and  Y(z). 

Let 

U(z)  - (Uj (z) , U2(z),  ...,  UR(z))  (51) 

where  U^z)  is  the  z-transform  of  u^(t)  which  is  the  upgoing  wave  in 
the  i-th  layer  of  the  complete  model.  Since  the  response  of  the  com- 
plete model  is  the  superposition  of  the  responses  of  primaries,  secon- 
daries, tertlaries,  etc.  models,  we  can  write 

00 

u(z)  - vz>  (52) 

where  U (z)  is  given  in  Eq.  (45b).  Notice  the  difference  between 
n 

U^tz)  and  IJ ^ ( z ) ; the  former  is  a scalar  function  which  is  associated 
with  the  i-th  layer  in  the  complete  model,  whereas  the  latter  is  a 
vector  which  denotes  all  the  upgoing  waves  in  the  i-aries  model  . 
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Substituting  Eq.  (41)  into  (52),  we  have 

ao 

U(z)  - £(HA2)n-1  U^z) 


n»  1 

where  the  explicit  dependence  of  H on  z is  omitted  for  notational 
simplicity.  Since  the  infinite  series  in  the  right  hand  side  of 

Eq.  (53)  converges  to  U(z),  we  can  write 
00 

.n-1 


Hence , 


or 


^.(HAj)"-1  - [l-HA2]-‘ 

U(z)  - [i-haJ-1  U^z) 


U^z)  - [l-HA2]  U(z) 


This  is  the  equation  for  the  upgoing  primary  waves  in  terms  of  the 
upgolng  waves  of  the  complete  model.  Let 

H - {h^}  i.j  - 1,2,. ...K  . 

Then  h^  is  the  corresponding  z-transform  of  h^  given  by  Eq.  (31); 
i.e. , 


hiJ  " hij 


2i  2 


(53) 


(54) 


(55) 


(56) 


(57a) 


(57b) 


Then,  substituting  (57  j)  into  the  vector  equation  (56)  and  taking  its 
first  component,  we  have 

Uifl(z)  - (l  + rQ  ^(z)  + rL  h12  U2(z)  + ...  + rR-1  h1R  UR(z) 


or 


Ui,i(z)  - Uj(z)  - [-rQ  fin  Uj(z)  - rl  hl2  U2(z)  - ...  - rk-1  \(z)] 

(58) 


Since  U.  , (z)  represents  the  primary  reflections  and  U, (z)  represents 
the  complete  response,  the  terms  in  the  brackets  in  Eq.  (58)  should 
represent  all  the  multiple  reflections.  Our  objective  is  to  suppress 
those  from  U^(z).  To  achieve  this  purpose  we  should  express  the 
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bracketed  terms  in  terms  of  ^(z)  and  our  known  quantities  U^(z),  r^ 
and  M(z)  (U^(z)  is  obtained  directly  from  Y(z)  by  Eq.  (3)).  But,  this 
is  impossible  because  those  terms  in  the  brackets  are  functions  of 
t^,t2,  • . . »T^  and  r^,r2,...,rK  which  we  do  not  know  and  moreover 
U2(z),  U^(z),  ...»  U^(z)  are  not  measurable.  However,  the  first  term 
in  the  brackets,  especially  h^,  which  is  also  a function  of  t^,t2,...,tk 
and  r^,r2>. . . ,rR,  can  be  expressed  in  terms  of  ^(z)  and  the  known 
quantities  from  Eq.  (46a),  and  this  fact  enables  us  to  make  an  approxi- 
mation of  the  primaries  U .(z)  from  U (z). 

1 » 1 1 

Taking  the  first  component  of  the  vector  equation  (46a),  we  have 


Hence , 


Ux  x(z)  - h^d  + r^  M(z) 


‘11 


Ul,l(2> 

(ld-r0>  M(z) 


(59) 

(60) 


Substituting  this  into  Eq.  (58), 

r ui  i(z)  1 

0.  . (z)  - U (z)  - [-r  ' --  U (z)  + a(z)  J 

i*1  1 L 0 (l+r0)  M(z)  1 J 


(61) 


where 


a(z)  - -rj  h12  U2(z)  - ...  - rK-1  h1R(z) 
From  Eq.  (61),  it  follows  that 


U1,1(Z) 


UL(z) 


a(z) 


_ rQ  Ul(2)  TJ_  ui(2) 

(1+rg)  M(z)  (1+rg)  M(z) 


From  Eqs.  (3)  and  (11),  we  have 

1 


U1(Z)  ' d-rj 


Y(z 


ui,i(z)  “ a^Yi(2)- 


(62) 


(63) 


(64) 

(65) 
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Substituting  these  into  Eq.  (63)  for  U.  (z)  and  U.  .(z),  we  get 

1 i » 1 

Y(z) 


Y^z) 


(l-rQ)  a (z) 


1 - 


Y(z) 


Y(z) 


(66) 


1 -■ 


(1-rJ)  M(z) 


(1-rJ)  M(z) 


Usually  the  second  term  in  the  right-hand  side  of  Eqs.  (66)  is  quite 

small  in  magnitude  compared  with  Y,(z),  and  especially  when  r >>  r 

1 o i 

(i  “ 1»2, . . . ,K) , as  in  most  geophysical  situations,  this  term  is  almost 
negligible.  In  this  case 


Let 


and 


Yj(z) 


Fj(z) 


8(z) 


Y (z) 


ro 

Y (z) 

d-ro) 

M(z) 

Y(z) 

ro 

Y(z) 

d-rj) 

M(z) 

-rQ)  a(z) 

ro 

Y(z) 

(i-rj) 

M(z) 

(67) 


(68) 


(69) 


Then  Eq.  (66)  becomes 

Yj  (z)  - Fl(z)  - 8(z) 
or 


Fj(z)  **  Yj  (z)  + 8(z) 


(70) 


Let  the  relation  given  by  Eq.  (68)  be  F^.  Then  F1  is  the  inverse 
filter  we  were  looking  for;  i.e.,  given  the  synthetic  seismogram  Y(z) 
with  knowledge  of  the  surface  reflection  coefficient,  r^,  and  the  input 
waveform,  M(z),  F^  suppresses  some  amount  of  (actually  the  most  signi- 
ficant portions  of)  the  multiple  reflections  from  Y(z)  and  yields  an 
approximation  of  the  primaries  Y^(z). 


If  the  input  is  an  impulse,  then  M(z)  ■ 1 and  Eq.  (68)  reduces  to 

Y(z) 

F.(z) (71) 

1 0 

1 V Y(Z) 

(1-rJ) 

In  this  case  Y(z)  represents  the  cransfer  function  of  the  layered  media. 

A simulated  result  is  shown  in  Figures  8 through  10  for  the  three 
layer  example  given  in  Section  III.  The  impulse  response  of  the  media  and  its 
filtered  output,  obtained  by  Eq.  (71) ? are  shown  together  in  Figure  8. 
Figure  9 depicts  the  input  waveform  and  Figure  10(a)  is  the  synthetic 
seismogram  due  to  this  input  convolved  with  the  impulse  response.  The 
result  of  applying  the  filter  in  Eq.  (68)  to  the  seismogram  is  shown 
in  Figure  10(b).  The  three  prominent  peaks  in  Figure  10(b)  represent 
the  primaries  Y^(z).  The  remaining  small  ripples  represent  the  additional 
term  B(z)  in  Eq.  (69)  or  (70).  We  see  that,  in  this  example,  significant 
multiple  reflections  are  almost  completely  eliminated  by  our  filter. 
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VII.  Another  Filter,  which  makes  use  of  Additional  Information  about 
r^  and  . 

Suppose  the  reflection  coefficient  r^  associated  with  the  1st  inter- 
face of  the  media  and  the  one  way  travel  time  of  the  1st  layer  are 
known  in  addition  to  knowledge  of  r^  and  the  input  waveform.  This  assump- 
tion may  be  applicable  in  the  marine  situation.  Making  use  of  this  extra 
Information  we  want  to  find  an  inverse  filter  which  suppresses  some 
additional  multiple  reflections. 

In  the  following  we  will  show  that  the  second  term  -r ^ h^  ^(z)  in 
the  brackets  in  Eq.  (58)  can  be  expressed  in  terms  of  U (z)  and  the 

A » 1 

known  quantities  rg»ri»T|  and  M(z). 

From  Eqs.  (31a),  (31b)  and  (57b),  we  have 


hij  qi~j 


for  i > j 


(72) 


Especially, 


or 


12 

'21 


‘12 


(1-rp 


qL  (1+rj) 


o-^)  . 

(I-HTj)  h21 


(73) 


(74) 


Now,  taking  the  second  component  of  the  vector  equation  (46a),  we 
have 
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P1.2<«> 

(l+r0)  M(z) 


(75) 


Hence , 


(1'rl)  U1.2(z) 
‘l2  " (l+r0)(l+rl)  M(z) 


(76) 
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Since  U.  »(z)  is  not  measurable,  we  still  need  an  expression  for  it. 

To  obtain  U.  _(z)  we  should  refer  to  the  state  equations  of  the  pri- 
l 

maries  model.  The  z-transforms  of  the  first  two  equations  of  Eq.  (10) 


are 


z nlD1  ^ (z ) - (1+Tq)  M(z) 


z"niUlfl(z)  “ rl  D1,1(Z)  + (1'rl)  ”l,2(z) 


(77a) 

(77b) 


From  these,  we  get 


-n, 


z U,  ,(z)  - (l+rQ)r1  z M(z) 


Ul,2(z> 


1,1 


(l-rx) 


(78) 


Substituting  Eq.  (78)  into  (76)  for  U,  .(z),  we  have 


z'ni  ci.i(z) 


ri z 


nl 


12  (l+rgJd+r!)  M(z)  (1+ri) 


(79) 


Rewriting  the  first  two  equation  in  Eq  (5)  in  z-transf ormation, 
-n 


z 1 (z)  » -rQ  (z)  + (l+rQ)  M(z) 

z 1 U1(z)  - rL  Dx(z)  + (1-r^  U2(z) 


(80a) 

(80b) 


From  these,  we  get 


U2(z) 


(z  1+rorl  z ^ Ul^  “ ^1+ro^rl  2 1 


(1-rj) 


(81) 


Thus,  the  first  two  terms  in  the  brackets  of  Eq.  (58)  are  expressible 

in  terms  of  known  and  measurable  quantities  and  U .(z).  We  write  Eq.  (58) 

1 , i 

as 


where 


U1  l^z)  * Ui<z)  “ [-r0  ^11  Ui^z^  ” ri  ^*12  U2(z^  + Y(z)]  , 


Y (z)  - -r2  h13  U„ (z)  - ...  - rK_L  h1R  UR(z) 


(82) 


(83) 
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Substituting  Eqs.  (60),  (79)  and  (81)  into  Eq.  (82),  we  have 


2 3 2n.  o 2n, 

(1~2rl-r0rlZ  > <l+r0)rj  z 1 

ui  ,(z)  * U,  (z)  + - — 1 M(z) 


1,1 


(1-r^) 


1 


(1-rJ) 


V!l] 

Ld-Tj)  (l+rQ)(l-rp  M(z)  J 1,1 


(z)  - y (z) 


(84) 


Finally,  solving  Eq.  (84)  for  ^(z),  we  have  the  following  equation  : 


U1,1(Z> 


,,  2 3 2n,  -i  2n, 

(l“2r , - rrtr , z l)  U,  (z)  + (l+rQ)r^  z 1 M(z) 


1 0 1 


r 


l - 


(rc+ri  2 ui(z) 


- <5  (z) 


(85) 


(l+r0) 


M(z) 


where  6(z)  is  a residual  term  which  is  deduced  from  Eq.  (84).  Substituting 

Eqs.  (64)  and  (65)  into  Eq.  (85),  we  obtain 

, 2 3 2n.  ? o 2n, 

(l-2r  -r  rlz  ‘)  Y(z)  + (l-r‘)r  z 1 M(z) 

Vz>* L—Q-i =2rT ^ n(z)  (86) 

(r  +r.  z L)  Y (z) 

1 - 

(1-rJ)  M(z) 


where 


n(z)  = ( l~r q ) <5 (z ) 


Now,  let 


F2(z) 


2 3 2n,  o o 2n i 

(1-2^-^  z L)  Y(z)  + (l-rg)rj  z 1 M(z) 


1 - 


(rQ+r1  z“2nl)  Y (z) 


(1-r2) 


M(z) 


(87) 


(88) 
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Then 


F2(z)  - Y^(z)  + n(z) 


Let  the  relation  in  Eq.  (88)  be  denoted  by  . Then,  F^  is  the 

inverse  filter  we  were  looking  for.  F2(z)  is  a better  approximation 

of  Y^(z)  than  F^(z)  given  in  Eq. (68).  For  the  simulation  purposes, 

our  expression  for  F2(z)  in  Eq.  (85)  is  not  so  useful  because  of  its 

complicated  structure.  In  the  following  section  we  will  show  that 

filter  F2  is  just  equivalent  to  the  successive  application  of  filter  F 

in  two  stages.  From  this,  a recursive  method  to  generate  F (z), 

n 

n-2,l,...,K,  which  are  the  successive  approximations  of  Y^(z)  with 
Fk(z)  = Y^(z),  will  be  developed  in  Section  IX. 


VIII.  Effects  of  Filters  and 


Let  R(z)  be  the  transfer  function  of  the  system  of  layered  media. 


It  is  given  by 


R(z)  ‘ 5(f> 


Then,  Eq.  (68)  for  filter  can  be  written  as 


Fl(z) 


1 V R(z) 

d-»S> 


If  the  input  is  an  impulse,  then  M(z)  - 1 and  Eq.  (91)  reduces  to 


F2(z) 


1 V R(z) 

d-rj) 


Now,  we  define  another  transfer  function  G(z)  as 

G(z)  - R(z)  + rQ  ( 

The  additional  term  r q represents  the  direct  reflection  of  the  impulse 
at  the  surface  at  time  zero  (Figure  11a).  Notice  that  R(z)  does  not 
include  this  direct  reflection  because  in  our  state  space  model  Y(z) 
was  defined  as  Y(z)*  (l-r^)  U^(z)  and  it  excluded  the  direct  reflec- 
tion term  Tq  M(z)  (see  Eq.  (3)).  The  expression  for  filter  F^  in 
terms  of  G (z)  is  obtained  just  by  substituting  Eq.  (93)  for  R(z) 
into  Eq.  (91)  or  (92).  In  a similar  manner,  we  define  R^(z)  and  G^(z) 
to  be  transfer  functions  of  the  subsystem  of  the  layered  media  below 


the  1st  layer  such  that 


G^  (z)  - Rx(z)  + rx 


This  is  depicted  in  Figure  lib.  Then,  observing  the  ray  diagram  in 
Figure  11c,  we  can  write 
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P(z) 


R(z)  _-nl 


(l-rQ) 


n, 


"1  r0  2 
Q(z)  =■  (1+rg)  z 1 - ^~r  j R(z) 


G (z)  » P(z) 
1U;  Q(z) 


From  these  we  get  the  relation  between  G^(z)  and  R(z), 


Gx(z) 


z 2ni  R(z) 


(l-rQ)  - rQ  R(z) 


or 


2n 


1 


. z * G (z) 
R(z)  - (1-0 


"(v  Zni 

l + r0z  1 Gl(z) 


Substituting  Eq.  (96b)  into  Eq.  (91),  we  have 

2 2n  l 

Fj(z)  - (l-rQ)z  1 Gj(z)  M(z) 

If  the  input  is  an  impulse,  then  M(z) = 1,  and 

2 2nl 

F^z)  * (l-rj)z  1 G^z) 


(95a) 


(95b) 


(95c) 


(96a) 


(96b) 


(97a) 


(97b) 


This  result  indicates  the  fact  that  the  output  of  the  filter 
F^(z)  is  the  same  as  the  output  of  the  system  of  the  layered  media  just 
below  the  1st  layer,  which  is  observed  at  the  surface.  This  is  a rather 
interesting  result  because  the  impulse  response  (or  transfer  function)  of  this 
subsystem  can  be  obtained  from  the  output  of  the  filter,  i.e.,  from  eg.  (97b) , 

_2nl 

Gl(z)  - -5 3-F  (*)  (98) 

(l-r0) 

This  result  also  indicates  that  filter  F^  eliminates  all  the 
multiple  reflections  which  are  ever  reflected  off  of  the  surface, 

When  the  reflection  coefficient  of  the  surface  is  relatively  large, 
the  multiple  reflections  removed  by  filter  F^  are  those  which  have 
significant  magnitudes  and  the  effect  of  the  filter  is  especially 
remarkable. 
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Filter  F 2 given  in  Eq.  (88)  is  expressed  in  terms  of  R(z)  as 

2 3 2n,  2 3 2n, 

(l-2r  -r  r,  z A)  R(z)  + (l-r^)r,  z 
F2(z) L_0_I 0_L M(z) 

(r0+rl  z > 

1 — S R(Z) 

(l-r$) 


(99) 


Substituting  Eq.  (96b)  for  R(z)  into  Eq.  (99),  we  get 

, 2n,  <l-2r?)  6 (z)  + r? 

F (z)  - (l-r‘)  z 1 1 M(z) 

° l-riG(z) 


(100) 


Again,  substituting  Eq.  (94)  into  Eq.  (100)  for  G^(z),  we  obtain  the 
following  equation. 


2 ^nl 
F2(z)  - (l-r£)  z 1 


rl  + 


RL(z) 


1 - 


(l-r‘) 


2 Rx(z) 


M(z) 


(101) 


Here,  R^(z)  is,  as  we  defined  in  Eq.  (94),  the  transfer  function  of  the 

subsystem  of  the  layered  media  below  the  1st  layer,  excluding  the 

direct  reflection  term  r^  at  time  zero.  The  second  term  in  the 

brackets  in  Eq.  (101)  is  exactly  of  the  same  form  as  filter  F^, 

hence  this  term  represents  the  filtered  output  of  the  subsystem 

2 2ni 

below  the  1st  layer.  The  first  term  (1-r^)  z r^  M(z)  in  Eq.  (101) 
represents  the  primary  reflection  associated  with  the  1st  interface 
of  the  media.  Hence,  F2(z)  can  be  obtained  by  applying  filter  F^ 
twice  successively,  first  to  the  given  system  and  then  to  the  sub- 
system below  the  1st  layer;  i.e., 

(1)  apply  F1  to  Y(z)  to  get  F^z), 

(2)  then,  compute  R^(z)  by  Eqs.  (98)  and  (94)  from  F^(z)  and  apply 
F^  again  to  R^(z),  F2(z)  is  then  obtained  by  substituting  this 
result  into  Eq.  (101). 

Now,  if  we  define  G2(z)  and  R2(z)  to  be  the  transfer  functions 
of  the  subsystem  below  the  second  layer  in  the  same  way  as  we  de- 
fined G^(z)  and  R^(z),  such  that 

G2(z)  - R2(z)  + r2  , (102) 
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then,  it  is  not  difficult  to  show  that 


R^z) 

V> 

(1-r^ 

Equation  (103)  is  analogous  to  Eq.  (97b).  Substituting  this  into  Eq. 

(101),  we  have 

r 2 2ni  2 2 2ni+2n?  _ 

F2(z)  - [(1-rJ)^  z + (l-rj)(l-rp  z G2(z)j  M(z)  (104) 

This  result  shows  that  F2(z)  is  just  the  output  of  the  subsystem  below 
the  2nd  layer,  which  is  observed  on  the  surface,  plus  the  primary  re- 
flections reflected  from  the  1st  interface.  This  is  described  in  Fi- 
gure 12.  This  result  also  indicates  that  filter  F^  eliminates  all  the 
multiple  reflections  which  are  ever  reflected  down  from  the  surface  and 
the  1st  interface. 


2 2n2 

(l-r‘)  z L G2(z) 


(103) 
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In  Section  VIII  we  have  already  recognized  that  F ^ (z)  can  be  obtained 
in  a recursive  way  applying  filter  in  two  stages,  first  to  the  whole 
systemR(z)  and  then  to  the  subsystem  R^(z)  which  is  obtained  directly  from 
the  result  of  the  previous  step.  This  procedure  can  be  extended  to  apply 
filter  F^  to  the  consecutive  subsystems  R 
F^ (z) ,F^ (z) , . . . ,FR(z) . The  subsystems  are  computed  recursively  during  this 
procedure.  The  outputs  F^ (z) .F^ (z) , . . . ,FK(z)  are  then  successive  approxi- 
mations to  the  primaries  Y^(z).  It  will  be  shown  that  in  the  deterministic 
situations  (i.e.,  perfect  measurements)  the  final  output,  FR(z),  represents 
the  pure  primaries  of  the  system. 

Observe,  that  to  compute  F^(z)  we  need  to  know  the  value  of  r^,  and 
to  obtain  F^Cz),  additional  knowledge  of  x^  and  r^  is  necessary.  Likewise, 
in  order  to  apply  the  recursive  procedure  to  obtain  the  third  output  F^(z), 
knowledge  of  x^  and  r^  is  required,  and  in  general,  to  get  Fn(z)  at  the 
n-th  stage  of  the  procedure  we  need  knowledge  of  x^  ^ and  rn  In  this 
paper  we  just  assume  that  those  values  are  known  in  each  step  of  the  pro- 
cedure (by  estimation  or  whatever).  In  a later  paper  we  shall  discuss 
the  estimation  of  these  quantities. 

Just  as  we  defined  G^(z),  R^(z)  and  G2(z),  R2(z)  in  Section  VIII,  let 
G^z),  R^z)  be  the  transfer  functions  of  the  subsystem  of  the  media  below 
the  i-th  layer;  G^(z)  Includes  the  direct  reflection  term  r at  time  zero 
and  Rj(z)  is  just  the  same  as  G^z)  except  that  it  excludes  the  direct  re- 
flection term;  i.e., 

Gt(z)  - Ri(z)  + rt  (105) 


2(z)  »R3(z) . • • • .^(Z)  generating 
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(107) 


Let  us  define  X^(z)  to  be 


X1(z)  - 


^(z) 


2TRi(2) 


(1-rp 


Then,  X^(z)  is  the  output  of  the  filter  which  is  applied  to  the 
subsystem  R^z).  Substituting  Eq.  (106)  into  (107)  for  R^z),  we  have 

Xt(z)  - (1-r*)  z2ni+1  G1+1(z) 

From  Eqs.  (105)  and  (108),  we  obtain 


(108) 


Ri+l(z)  ” ,,  2.  Xi(z)  " rd 


(109) 


(1-rp 


This  equation  states  the  important  fact  that  subsystem  Ri+^(z)  can  be 

obtained  from  the  filtered  output  of  subsystem  R^(z).  This  Eqs.  (107) 

and  (109)  represent  a recursive  relation  for  Ri(z)'s,  i-l,2,...,K, 

with  R(z)  as  the  starting  value.  Equations  (97a)  and  (104)  can  be  written 

in  terms  of  R^(z)  and  R2(z)  as 

[2  2ni  2 ^nl  T 

(l-r^rj  z 1 + (1-rJ  z 1 Rj(z)|m(z)  (110) 

and  „ 

[2  zni  2 2 2n]+2no 

(l-rQ)r1  z + (1-r^) (l-r1)r2  z 

2 2 2n.+2n 

+ (l-rj)(l-rj)  z R2(z)  I M(z)  (111) 


Observing  the  expressions  for  F^(z)  in  Eq.  (101),  we  can  see  that 
F2(z)  is  obtained  if  we  replace  R^(z)  in  Eq.  (110)  by  X^(z),  the 
filtered  output  of  R^z).  Likewise,  we  expect  that  F^(z)  will  be 
obtained  if  we  replace  R2(z)  in  Eq.  (Ill)  by  X2(z).  This  proce- 
dure can  be  continued  to  generate  F^ (z) ,F^ (z) , . . . ,F^(z)  with  the 
help  of  the  recursive  relations  given  by  Eqs.  (107)  and  (109). 
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Then,  Fi+j(2)  has  the  following  structure; 


W'> 


r i / i- 

E(n 

_i-l  \ k-l 


, 2n,  , 

n d-r?)  z k+1 
k-0 


)r1+(n 

• Vk-o 


d-r^)  z 


))xi<2>] 


(112) 


The  first  term  in  the  brackets  represents  the  primaries  which  are 

associated  with  the  first  n interfaces  (interfaces  Hi)  of  the  media 

and  the  second  term  represents  the  filtered  output  X (z)  of  the  sub- 

n 

system  Rn(z)  observed  at  the  surface.  This  is  depicted  in  Figure  13. 
Since  the  effect  of  the  filter  is  to  remove  all  the  reflections  which 
are  ever  reflected  down  off  the  top  surface  of  the  media  (subsystem), 
Fn(z),  as  the  result  of  the  n-th  stage  of  the  procedure  is  free  from 
all  those  multiple  reflections  which  are  ever  reflected  down  off  the 
surface  and  the  first  (n-1)  interfaces.  If  we  continue  this  procedure, 
after  the  K-th  stage  all  the  multiple  reflections  will  be  removed 
completely  from  the  seismogram  and  only  primary  reflections  will  remain. 
To  see  this,  observe  from  Eq.  (112)  that 


fr(z) 


EG  I 


2.  2”k+l 

«-rk)  z 


H-iw] 


But  from  Eq.  (108) 


XK-l(z)  “ *1_rK-P  2 GK(z) 


(113) 


(114) 


where  the  subsystem  G^(z)  below  the  K-th  layer  is  just  the  basement 


gk(z)  - rk  . 


Hence, 


WsG (l<)  *2”k+1H 


(115) 


(116) 


which  represents  the  pure  primaries  of  the  system. 
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Summarizing  the  recursive  procedure  : 
Starting  Equation 


Vz) 


R(z) 


1 V R(z) 


(1-r^) 


Recursive  Relations 


“2n . 


Rt(z) 


(1-r 


2 Xi-l(z)_ri 


i+1) 


Xi(z) 


Rt(z) 


1 TRi(2) 

(1-^) 


(117) 


(118) 


(119) 


Output  Equations 


Fx(z)  - XQ(z)  M(z) 


(120) 


(1-r^) 


1 2 
n (1-r") 

k-0 


where  i ■ 1,2, . . . ,K-1. 


(121) 


Figure  14  depicts  filtered  outputs  F^(z)  and  F^(z)  for  the  three 
layer  example  which  are  obtained  by  the  recursive  procedure  above 
(This  is  continuation  of  Figure  10).  The  last  result,  F^(z),  which  is 
shown  in  Figure  14(b)  exhibits  the  three  pure  primaries  of  the  system. 
Since,  in  this  example,  the  considerable  amount  of  multiple  reflections 
was  removed  already  in  the  first  stage  (see  Figure  10(b)),  F^(z)  and 
F^(z)  in  Figure  14  shows  little  changes  in  improvement. 


34 


This  is  due  to  the  fact  that  the  surface  reflection  coefficient,  r^, 
in  this  example  is  much  greater  than  r^,r2  and  r^  in  magnitudes.  Fi- 
gure 15  depicts  another  example  of  the  recursive  filtering.  Figure 
15(a)  is  the  complete  response  of  a six  layer  media.  The  values  of 
the  reflection  coefficients  and  the  one-way  travel  times  of  the  media 
used  in  this  example  are  ; 


o 

u 

- 0.68  , 

rl  - 

0.40 

’ r2 

—0.32  , 

r3  - 0.53, 

r4 

—0.78  , 

r5  * 

0.71 

’ r6 

- 0.65  , 

and 

T1 

- 0.07  , 

T2  " 

0.04 

’ x3 

* 0.115, 

- 0.09, 

t5 

- 0.035, 

t6  * 

0.13 

9 

The  same  input  given  in  Figure  9 is  used  in  this  example.  Figures  15 
(b)-(g)  show  successive  filtered  outputs,  F^z),  F2(z),  ....  F&(z), 
which  are  generated  by  our  recursive  filtering  procedure.  Again, 

F^(z)  in  Figure  15(g)  represents  the  primaries  of  the  six  layer  media. 
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X.  Conclusions 

We  have  presented  a simple  Inverse  filter  which  suppresses  a fair 
amount  of  the  multiple  reflections  in  a synthetic  seismogram.  This  filter 
requires  knowledge  of  the  surface  reflection  coefficient  and  the  input  wave- 
form. The  Bremmer  series  decomposition  played  a key  role  in  its  develop- 
ment. The  filter  was  shown  to  be  especially  useful  when  the  surface  re- 
flection coefficient  is  relatively  large  (as  in  most  geophysical  situa- 
tions) in  which  case  significant  multiple  reflections  are  almost  completely 
removed  so  that  the  output  is  a good  approximation  of  the  primaries  of  the 
layered  system. 

The  recursive  filtering  method  presented  herein  demonstrates  the 
possibilities  that  the  subsystems  of  the  layered  media  can  be  revealed 
from  the  seismogram  by  applying  the  filter  successively  and  that  pure 
primaries  of  the  system  can  be  obtained  thereof.  The  actual  application 
of  this  recursive  procedure  requires  the  estimation  of  the  reflection 
coefficients  and  the  oneway  travel  time  in  each  stage  to  perform  the  next 
recursion;  but,  this  in  turn,  provides  motivation  to  use  this  procedure 
in  estimating  those  quantities.  The  estimation  scheme  is  now  under 
research.  Also,  work  is  in  progress  on  developing  a version  of  this 
filter  that  is  applicable  to  noisy  seismogram  or  realistic  field  seismic 
data  and  to  the  situation  when  the  source  waveform  is  not  available. 
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